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Abstract 

Multiscale modelling methodologies build macroscale models of 
materials with complicated fine microscale structure. We propose a 
methodology to derive boundary conditions for the macroscale model 
of a prototypical non-linear heat exchanger. The derived macroscale 
boundary conditions improve the accuracy of macroscale model. We 
verify the new boundary conditions by numerical methods. The tech¬ 
niques developed here can be adapted to a wide range of multiscale 
reaction-diffusion-advection systems. 
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1 Introduction 

Multiscale modelling techniques are a developing area of research in engi¬ 
neering and physical sciences. These techniques are needed when the system 
being modelled possesses very different space-time scales and it is infeasible 
to simulate the whole domain on a microscale mesh (Dolbow et al. 2004, 
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Figure 1: A schematic diagram of a heat exchanger. The red pipeline carries 
fluid to the right, and the blue pipeline carries fluid to the left. Heat exchanges 
between the pipes. 


Kevrekidis & Samaey 2009, Bunder & Roberts 2012). Macroscale boundary- 
conditions are rarely derived systematically in such works. Instead macroscale 
boundary condition are often proposed heuristically (Pavliotis & Stuart 2008, 
Mei & Vernescu 2010, Mseis 2010). We developed a systematic method 
to derive boundary condition for one-dimensional linear problems with hne 
structure by cell mapping (Chen et ah 2014). Here we extend the method to 
a prototypical non-linear heat exchanger problem. 

We mathematically model the counter flow two-stream heat transfer shown 
in Figure 1. Let x measure nondimensional distance along the heat exchanger 
which is of length L, 0 < x < L, and let t denote nondimensional time. The 
held a(x, t) is the temperature of the huid in one pipe and held b{x, t) is that 
in the other pipe. A quadratic reaction is included as an example nonlinearity 
to give the nondimensional microscale pdes 
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The lateral diffusion also included in these pdes makes the derivation of 
macroscale boundary conditions challenging. 

Various mathematical methodologies derive from pdes (1) the macroscale 
model 

f = + + + (2) 

for the mean temperature C{x, t) := [a(a;, t) + b{x, t)]/2 . For example, centre 
manifold theory rigorously derives this effective macroscale model (Roberts 
2013), as does homogenization (Pavliotis & Stuart 2008, Mei & Vernescu 
2010). This macroscale model combines an effective cubic reaction, with an 
effective nonlinear advection, and enhanced lateral diffusion. The necessary 
analysis to derive the model (2) is based around the equilibrium a = b = 0, 
and applies to the slowly-varying in space solutions in the interior of the 
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domain. For example, in the interior it predicts the temperature helds are 



+ OiC^ + dl). 


(3) 


The challenge of this article is to provide sound boundary conditions for the 
macroscale model (2). 

Prototypical microscale boundary conditions for microscale system (1) are 
taken to be the Dirichlet boundary conditions 


a(0,f) = ao, b{0,t)=bo, a{L,t) = a^, and b{L,t) = bL, (4) 


where oq, ul, bo and b^ are potentially slowly varying functions of time. 
Section 3 derives nonlinear macroscale boundary conditions (12) for the 
macroscale mean temperature model (2). For example, the linearisation of 
the macroscale boundary condition (14) derived from the Dirichlet boundary 
conditions (4) is the Robin condition 
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at X = 0 . 


Importantly, this is not a Dirichlet boundary condition despite the microscale 
boundary conditions and the dehnition C := {a + b)/2 together suggesting 
boundary conditions are (incorrectly) C{0,t) = (oq + bo)/2. 

Figure 2 plots microscale and macroscale solutions for the heat exchanger 
at a particular time. The two solid lines plot the microscale solution a{x,t) 
and 6(x, t) of microscale pde (1) with microscale boundary conditions (4). The 
black dashed line plots the mean temperature model (2) with classic Dirichlet 
boundary conditions C{0,t) = {ao + bo)/2 and C{L,t) = (oi, + 6 l)/ 2 as 
would be commonly invoked (Mei & Vernescu 2010, Mseis 2010, Ray et ah 
2012). The macroscale model (2) performs poorly with these heuristic Dirichlet 
boundary conditions, especially in the interior of the domain (here 5 < x < 25). 
But the interior is where the macroscale model (2) should be valid. The 
macroscale model (2) represents the interior dynamics but cannot resolve 
the details of boundary layers (Roberts 1992). With our derived boundary 
conditions, the macroscale solution (red line in Figure 2) hts the microscale 
solution (solid lines) in the interior: the microscale helds a{x, t) and b{x, t) 
being given by equation (3). Our systematic derivation of boundary conditions 
is needed for macroscale models to correctly predict the interior dynamics. 

The key to our approach is to explore the effect of boundary layers by 
treating space as a time-like variable (Chen et ah 2014, e.g.). However, the 
heat exchanger problem (1) is challenging because of the nonlinearity. Here, 
a normal form coordinate transformation separates the spatial evolution 
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Figure 2: Example solutions of the heat exchanger (1) in domain 0 < x < L = 
30 at time t = 21. The two solid lines plot the the temperature of the two pipes, 
a{x,t) and b{x,t). The dashed lines are solution of the macroscale model (2) 
at t = 21: (black dashed) with heuristic Dirichlet boundary conditions; and 
(red dash-dots) with our systematically derived boundary conditions. 


in the boundary layers into a slow manifold, stable manifold and unstable 
manifolds. This separation empowers a transformation of the given physical 
boundary conditions (4) into boundary conditions (13) for the macroscale 
interior model (2). 

2 A normal form of the spatial evolution 

The macroscale model (2) is slow so the dominant terms in the boundary layers 
are due to the derivatives of spatial structure. Thus, to derive macroscale 
boundary conditions for slow evolution (2) we treat the time derivative d/dt 
as a negligible operator (Roberts 1992). To put heat exchanger system (1) 
into the form of a dynamical system in time-like variable x we define a' 
and b' := Then rearranging system (1) in dynamical system form, with 
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dt = 0 for quasi-steady solution, gives 
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(5) 


We analyse the (spatial) dynamics of this system with ‘initial condition’ at 
X = 0 of the given microscale Dirichlet boundary conditions (4). 

Start by basing the analysis of (5) around the equilibrium at the origin, 
a = h = a' = b' = 0. The eigenvalues of the system linearised about the origin 
are 0 (twice) and ±-\/2 • The eigenvalues of zero corresponds to an eigenvector 
of (1,1,0, 0) and a generalised eigenvector of (—1,1,1,1). Hence the spatial 
ODE system (5) contains two centre (slow) modes, one stable mode and one 
unstable mode. 

Roberts (2014a, &) provides a web service to construct by computer algebra 
a coordinate transform which separates stable, unstable and centre manifolds. 
However, the web service does not directly apply to systems whose linearisation 
has a generalised eigenvector. To circumvent the generalised eigenvector, we 
choose to embed the ODE system (5) as the e = 1 member of the one parameter 
family of systems 
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where the last vector on the RHS is treated as a perturbative term: the 

parameter e counts the order of artihcial linear perturbation. The linear 

operator in system ( 6 ) now has no generalised eigenvector: its eigenvalnes 
are 0 (twice) and ±|, with corresponding eigenvectors (1,1, 0, 0), (—1,1,1,1), 
(—1,|,0, l), (—1,|,1,0). The web service (Roberts 2014a) then hnds a 
normal form coordinate transform as a multivariate power series in variables Sj 
and parameter e. Snbstituting e = 1 into the results reveals the centre 
manifold, stable manifold and unstable manifold for the spatial ODE (5). 

The three manifolds can be parametrised as we choose. We choose the deh- 
nition of the two parameters for the slow manifold to be the mean temperature 
Si = C and its spatial derivative S 2 = ^ ■ 

5 3 := I (3a — 36 — 3a' -|- 96'), 

5 4 := I (3a - 36 + 9a' - 36'), (7) 


51 :— ^ (a -|- 6 ) — C*, 

52 :- 2 (a +6) - —, 
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where S 3 parametrise the stable manifold, and S 4 parametrise the u n stable 
manifold. Then the web service (Roberts 2014a) derives the coordinate 
transform ( 8 ) giving a, b, a' and b' as a power series of Si, S 2 , S 3 , S 4 , and e: 

CL ~ Si — 'S2 T O.25S3 T 1.5s^ T 6S2 — ITsi-^s — 3.4s2S3 — O.OSSsg 


+ 0 . 75 s 4 + O.74S2S4 + O.56S3S4 — O.25S4 , ( 8 a) 

b ^ Si + S2 — 0.75s 3 — 1.5si — 6S2 — O.74S2S3 + O.25S3 

— 0.25 s 4 — I.IS1S4 + 3.4S2S4 — O.56S3S4 — O.O35S4 , (8b) 

a' S 2 + 1.5siS2 — 0.17s3 + 0.56siS3 + O. 9 IS 2 S 3 + O. 47 S 3 

+ 0.5s4 — 0 . 56 siS 4 + I.2S2S4 — O.33S4 , ( 8 c) 

b ~ S2 — 1.5siS2 T 0.5s3 T 0.56siS3 T I.2S2S3 — 0.33s3 

— 0.17 s 4 — 0.56 siS 4 + O.9IS2S4 — O.O47S4 , (8d) 


For simplicity we only record these and later expressions correct to quadratic 
terms in sj, that is, with cubic errors in the multinomial, and for simplicity we 
record coefficients to two significant figures, and here we evaluate the power 
series at e = 1 to recover a coordinate transform applicable to the original 
spatial system (5). The corresponding evolution of the spatial system (5) in 
these new variables Sj is also provided by the web service which determines 
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S2 , 

(9a) 

1.55i52, 

(9b) 

—0.6753 ~ 0.755351 — 0.945352 , 

(9c) 

TO .6754 — 0.755454 T 0.945452 . 

(9d) 


The normal form of the transformed system (9) has useful properties. Since 
^ = 5 ' 3 (si, 52)53 and ^ = ^' 4 ( 51 , 52)54 for some functions gj indicates that 
three invariant manifolds of the system (9) are 53 = 0, 54 = 0 and 53 = 54 = 0. 
From the linearisation of (9) these are the centre-unstable, centre-stable, and 
slow manifolds respectively. Further, because ^ and ^ are functions of 
only 5i and 52 , the planes of 5i and 52 constant are isochrons of the slow 
manifold (Roberts 1989) (sometimes called the leaves of the foliation, hbres, 
a hbration, hbre maps or fibre bundles (Murdock 2003, pp.300-2, e.g.)). 

One might query whether the transformation ( 8 ) and (9) is valid given 
that it is obtained by a power series in artihcial parameter e that is then 
evaluated at e = 1. The coefficients appear to converge well to the given 
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values, but as an independent check we also embedded the spatial ode (5) 
into the different family of problems 
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( 10 ) 


Performing the same algebraic construction, but from this quite different base, 
we hnd system (10) results in the same transform (8) and evolution (9). This 
confirms the perturbative approach via embedding. 


3 Projection reveals boundary conditions 

This section focuses on the boundary layer near x = 0. As shown by the 
solid lines in Figure 2, the microscale boundary conditions at x = 0 force a 
boundary layer in the microscale model (1). However, the macroscale model 
(2) does not resolve the boundary layer. Forcing the macroscale model to 
pass through (oq + bo)/2 introduces an error in the interior of the domain, 
as shown by the dashed blue line in Figure 2. Here we derive an improved 
boundary condition at x = 0 which reduces the interior error caused by poorly 
chosen macroscale boundary condition. 

The boundary layer must lie in the centre-stable manifold S 4 = 0 because 
if there was any component S 4 7 ^ 0 then this would grow exponentially quickly 
in space and dominate the solution across the whole domain. Algebraically we 
obtain the centre-stable manifold by substituting S 4 = 0 into the coordinate 
transform (8): the terms in (8) are arranged so that this simply means 
omitting the second line of each of the four pairs of lines. 

Then, as plotted schematically in Figure 3, the two Dirichlet boundary 
conditions (4) at x = 0 form a one dimensional curve (solid blue line) of 
allowed values in the three-dimensional centre-stable manifold parametrised 
by Si, S 2 and S 3 . Recall Oq and Bq are the boundary values at x = 0 from 
boundary conditions (4). The first two components on the centre-stable 
manifold (S4 = 0) of (8) reveal the microscale constraints on the boundary, 
upon defining s° := for i = 1, 2, 3 , 

s? - s^ + 0.25s|] + l.Ssf + 6 sf - l.ls?s[] - 3.4s^s^ - 0.035sf' 
s? + s^ - 0.75s[| - 1.5sf - 6 sf - 0.74s^s|] + 0.25sf 

( 11 ) 

These equations implicitly determines the solid blue curve in Figure 3. To 
explicitly describe the curve, recall that this is a power series with cubic 
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Figure 3: schematic plot of centre-stable manifold near the boundary a; = 0. 
The green plane is the centre manifold. The blue solid line is the set of values 
allowed by the microscale boundary condition at a; = 0. The blue dotted line 
is the projection of the microscale boundary values onto the slow manifold. 


errors and so we just need to consistently revert the series to give, say, the 
boundary values s? and Sg as a function of 5°, oo and Bq. Algebra determines 


s? ^ (0.256o-0.296g + 0.75ao-0.63ao6o + 0.18a^) 


+ s^ (0.5 - 2.860 + 3.7ao) + 3sf, 

( 12 a) 

Sg (—60 — 0.196o + flo — 2.3ao6o — 0.56ao) 


+ s^ (2 - 4.660 + 3.8ao) - 5.2sf. 

( 12 b) 


Since the slow dynamics in the interior of the domain must lie on the slow 
manifold Sg = 0 , appropriate boundary conditions for the interior dynamics 
must come from projecting these allowed boundary values onto the slow 
manifold. Because of the special normal form of the transformed system (9), 
the slow variables si and S 2 evolve independently of the fast variables S 3 
and S 4 , and the appropriate projection is the orthogonal projection along the 
isochrons Sg and S 2 constant onto the plane Sg = 0 —shown by the red lines 
in Figure 3. Equation (12a) describes the projected curve in the SiS 2 -plane 
illustrated by the blue dashed line in Figure 3. Recall from the amplitude 
dehnition (7) that C and si are the same. Hence substituting s? = C and 
S 2 = ^ into equation ( 12 a) forms the boundary condition at a; = 0 

dC 

C - (0.5 - 2.860 + 3.7ao) -3 

ox 

(0.256o - 0.296^ + 0.75ao - 0.63ao6o + 0.18ag) . 





(13) 






















4 A numerical example 


9 


This nonlinear Robin bonndary condition prodnces the correct macroscale 
slowly varying interior domain solutions of the microscale model pde (1). 

4 A numerical example 

As an example, let the boundary values be oq = 0.2f(t) and bo = 0 for 
f{t) = tanh^t varying smoothly but quickly from /(O) = 0 to 1. Macroscale 
boundary condition (13) gives the macroscale boundary condition at a; = 0 
for mean temperature model (2) 

dC fdC\^ 

C - [0.75/ + 0.5] — - 3 f — J = 0.15/ + 0.007^. (14) 

Macroscale boundary conditions on the right One method to derive 
the macroscale boundary conditions at a: = L is to appeal to symmetry. 
Dehne a new spatial coordinate x = L — x measuring distance from the 
boundary into the interior, and dehne new held variables d{x,t) = —b{x,t), 
b{x,t) = —a{x,t) and therefore C{x,t) = —C{x,t). Then the pde system (1) 
is symbolically identical in the tilde and plain variables. But the boundary 
conditions (4) at the right-boundary x = L are transformed to Dirichlet 
boundary conditions at x = 0 of a(0, f) = —bi and b{0,t) = —ai. Then the 
derivation of Sections 2 and 3 apply in the same way to the tilde problem. 
After computing the macroscale boundary conditions in coordinate x we 
transform back to the original coordinate x. 

For example, assume = 0 and bi = 0.2. The iteration scheme in 
Section 3 computes macroscale boundary condition on the boundary x = L 
(x = 0) 


-C- [0.75/ - 0.5] 


dC 



0.15/- 0.007// 


(15) 


By the chain rule || = —1, and substitute C{x,t) = —C{x,t) into boundary 
condition (15) 


C - [0.75/ - 0.5] 


dC 

dx 


+ 3 



0.15/- 0.007// 


(16) 


Numerics verifies the macroscale boundary conditions derivation 

Figure 2 plots a snapshot of the simulations on microscale model (1) and mean 
temperature model (2) for two cases: the Dirichlet boundary conditions Co = 
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(ao + bo) /2 and Cl [cll + bi) /2; and our systematic boundary conditions (14) 
and (16). Using finite differences we convert the system of two pdes (1) into 
a system of odes. Then Matlab’s odel5s applies a variable order method to 
compute the solution of the system of odes (Shampine et ah 1999). 

The numerical result is as expected. The macroscale model with systematic 
boundary conditions (13) model the interior domain microscale dynamics 
much better than that with heuristic Dirichlet boundary conditions. 

5 Conclusion 

We systematically derived macroscale boundary conditions from microscale 
Dirichlet boundary conditions. This methodology can be extended to mi¬ 
croscale Neumann and Robin boundary conditions. For the microscale 
Dirichlet boundary conditions, we evaluated the first two components of 
the centre-stable manifold (8a)-(8b) at a; = 0 to reveal the microscale bound¬ 
ary constraints (11). If the microscale boundary conditions were Neumann, 
we would use the last two components, (8c)-(8d). If the microscale boundary 
conditions were Robin, we would use linear combinations of the transform (8). 
The methodology also applies to more general multiscale modelling of pdes 
(Roberts 1992). 

Acknowledgements CC thanks Dr. Tony Miller for his advice and useful 
discussion, and CSIRO for their support in funding to participate in conferences 
and workshops. 
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